#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>

void GaussLegendre();
int frequency_of_primes (int n);

int main()
{   
   GaussLegendre();
  
   return 0;
}

void GaussLegendre(){
  
  FILE *fp = fopen("gausslegendre_sequencial.txt", "wt");
  int i;

  mpf_t a0,b0,t0,p0;
  mpf_t a1,b1,t1,p1;
  mpf_t pi,aux,aux1;

  if (fp == NULL) // Se nào conseguiu criar
  {
     printf("Problemas na CRIACAO do arquivo\n");
     return;
  }

  //definindo a precisao
  mpf_set_default_prec(10000000);

  mpf_init(a1);
  mpf_init(b1);
  mpf_init(t1);
  mpf_init(p1);

  mpf_init(pi);
  mpf_init(aux);
  mpf_init(aux1);

  // a0 = 1;
  mpf_init_set_ui(a0,1); 

  // b0 = 1/sqrt(2);
  mpf_init_set_ui(b0,2);
  mpf_sqrt(b0,b0); 
  mpf_ui_div(b0,1,b0);

  // t0 = 0.25;
  mpf_init_set_d(t0,0.25);

  //p0 = 1;
  mpf_init_set_ui(p0,1);

  
  for(i=0;i<30;i++){
     
     //a1 = (a0 + b0)/2;
     mpf_add(aux,a0,b0);
     mpf_div_ui(a1,aux,2);    

     //b1 = sqrt(a0 * b0);
     mpf_mul(aux,a0,b0);
     mpf_sqrt(b1,aux);

     //t1 = t0 - (p0*(pow((a0 - a1),2)));
     mpf_sub(aux,a0,a1);
     mpf_pow_ui(aux1,aux,2);
     mpf_mul(aux,p0,aux1);
     mpf_sub(t1,t0,aux);

     //p1 = 2*p0;
     mpf_mul_ui(p1,p0,2);

     //pi = pow((a1+b1),2)/(4*t1);
     mpf_add(aux,a1,b1);
     mpf_pow_ui(aux1,aux,2);
     mpf_mul_ui(aux,t1,4);
     mpf_div(pi,aux1,aux);
     
     // atualizando variaveis
     mpf_set(a0,a1); //a0 = a1;
     mpf_set(b0,b1); //b0 = b1;
     mpf_set(t0,t1); //t0 = t1;
     mpf_set(p0,p1); //p0 = p1;

  }

   mpf_out_str(fp,10,10000000,pi);
   fclose(fp);

   mpf_clear(a0);
   mpf_clear(b0);
   mpf_clear(t0);
   mpf_clear(p0);

   mpf_clear(a1);
   mpf_clear(b1);
   mpf_clear(t1);
   mpf_clear(p1);

   mpf_clear(pi);
   mpf_clear(aux);
   mpf_clear(aux1);

}
